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The set of coupled equations for the self-consistent propagator and the field expectation value is 
solved numerically with high accuracy in Euclidean space at zero temperature and in the broken 
symmetry phase of the 4> 4 model. Explicitly finite equations are derived with the adaptation of the 
renormalization method of van Hees and Knoll [12] to the case of nonvanishing field expectation 
• value. The set of renormalization conditions used in this method leads to the same set of countert- 

' erms obtained recently by Patkos and Szep in [21]. This makes possible the direct comparison of the 

accurate solution of explicitly finite equations with the solution of renormalized equations containing 
counterterms. The numerically efficient way of solving iteratively these latter equations is obtained 
by deriving at each order of the iteration new counterterms which evolve during the iteration process 
towards the counterterms determined based on the asymptotic behavior of the converged propaga- 
tor. As shown at different values of the coupling, the use of these evolving counterterms accelerates 
the convergence of the solution of the equations. 
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^ r-| I. INTRODUCTION AND MOTIVATION 

' The ^-derivable approximation [1, 2] (also called the two-particle-irreducible (2PI) approximation) represents a 
^ , self-consistent and systematically improvable approximation scheme, successfully applied in both out-of-equilibrium 
and equilibrium contexts. In the first case, it overcomes the secularity problem encountered with methods not using a 
self-consistent propagator (see the references given in [3]) and makes possible the description of late time dynamics of 
scalar [4-7] and fermionic [8] quantum fields. In the second case, it represents an adequate framework for the study of 
the phase transition in quantum field theories, where some systematic partial resummation of the perturbative series 
has to be implemented. This is needed because the perturbative expansion around the free theory cannot capture 
the development of collective phenomena and the change in the physical spectrum occurring at high temperature, 
where the temperature could compensate the smallness of the coupling [9] . More sophisticated resummations can be 
. . ' obtained by combining the 2PI formalism with an expansion in number of flavors [10, 11]. 

In the past few years much attention was devoted to understand the renormalization in the 2PI formalism. In [12] 
the temperature dependent subdivergences of the self-consistent propagator equation were localized in the symmetric 
phase of the 4> A model by expanding the propagator and the self-energy around the corresponding zero temperature 
quantities. It was shown that the divergences are resummed by a Bethe-Salpeter-type equation and that the self- 
energy can be made finite if one renormalizes both a four-point function satisfying the Bethe-Salpeter equation and its 
kernel. Imposing renormalization conditions, the renormalization of the self-consistent propagator equation at C(A 2 ) 
truncation of the 2PI effective action of the above model was performed in the symmetric phase in [13]. In [12] no 
attention was paid to subdivergences appearing at zero temperature, although as shown in [14, 15], the same structure 
of subdivergences pointed out in [12] appears already there. The use of a finite number of counterterms in [14, 15] 
made transparent that once the theory is renormalized at T = it remains finite at finite temperature, as expected. 
The extension of the renormalization method of [12] to theories with spontaneously broken symmetries was treated 
in [16], where it was explicitly applied to the O(N) model at Hartree level truncation of the 2PI effective action. 
In [17] it was shown for a generic truncation of the 2PI effective potential that the renormalization at nonvanishing 
field expectation value can be achieved by renormalizing the theory in the symmetric phase. The quantities which 
need to be renormalized in the symmetric phase to render the propagator and field equations finite are the various 
two- and four-point functions and the 2PI kernels of the Bethe-Salpeter equations which are satisfied by the former 
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quantities and which resum different types of subdivergences. These quantities arise naturally by taking sequential 
derivatives of the so-called 2PI-resummed effective action with respect to the field. The removal of the divergences by 
appropriate counterterms is achieved by imposing renormalization and consistency conditions 1 on the 2PI kernels and 
the two- and four-point functions. This diagrammatic renormalization procedure was extended to models containing 
fermionic [19] and gauge degrees of freedom [20], and the applicability of the 2PI renormalization method in the study 
of the time evolution of out-of-equilibrium scalar fields was shown in [7] . Working in the broken symmetry phase of 
one- and multicomponent scalar models, an alternative method for obtaining the counterterms up to skeleton order 
C(A 2 ) truncation of the 2PI effective action was given in [21, 22]. The renormalization of the O(N) model in the 1/N 
expansion of the 2PI effective action was achieved in [3, 17, 23] using the original fields and also the auxiliary field 
method. A renormalization method was developed in [24, 25] which, instead of considering the divergences of the 
usual Feynman diagrams and the Bethe-Salpeter equations resumming them, uses a specific resummation scheme at a 
given temperature where the model is parametrized. At a different temperature the method uses matching conditions 
at asymptotic momenta. 

In an equilibrium setting and beyond the Hartree approximation there are relatively few papers reporting on 
numerical solutions of renormalized 2PI equations, even in scalar models. In some phenomenological studies [26, 27] 
the renormalization is not done, e.g. in [26] the O(N) model is solved in the Minkowski space with some approximations 
by taking into account momentum-dependent corrections only in the imaginary part of the self-energy and treating 
the tadpole with a UV cutoff. 

The application of the 2PI method did not went yet beyond the demonstration of its features in the simplest models 
(e.g. the TV-component real scalar model with O(N) symmetry, in many cases restricted to N = 1). The explicitly 
finite self-consistent propagator equation of the </> 4 model was solved including the setting-sun diagram at vanishing 
field expectation value and at finite temperature in [13]. The renormalized O(N) model was solved at zero temperature 
and at vanishing field expectation value within the bare- vertex approximation of the auxiliary field formalism in [3] . 
The pressure of the one-component <j> model was calculated in Euclidean space and finite temperature in [28] solving 
a renormalized 3-loop effective action. Solving the renormalized field and propagator equation of this model in 
Minkowski space the phase transition was studied in [29]. It was found that including the field-dependent two-loop 
skeleton diagram at the level of the 2PI effective potential results in a second order phase transition [29]. 2 The study 
of thermal properties of the spectral function of the renormalized 4 theory in the symmetric phase was reported in 
[25]. 

However, especially in the above listed finite temperature studies, little technical details are given on the numerics, 
in particular, the rate of convergence of the iterative steps leading to self-consistent solutions. This fact does not allow 
us to easily infer the accuracy of the methods used in obtaining the results. In our opinion, the lack of standardized, 
well-tested algorithms might be the main reason for the rather scarce number of applications of the 2PI-approximation 
in the study of the thermodynamics of relevant field theoretical models. Our aim in the present work is to quantify the 
accuracy of iterative solutions and to compare different algorithms. We also would like to make connections between 
different renormalization methods. Our highly accurate solutions could serve as benchmark for more powerful methods 
to be used for a finite temperature solution of these equations. 

The outline of the paper is as follows. In Sec. II we introduce the model and derive the coupled set of self- 
consistent gap and field equations. In Sec. Ill we discuss their renormalization by adapting to broken symmetry 
phase the renormalization method of [12, 16]. We isolate the quantities which need to be renormalized and imposing 
renormalization conditions we recover, on the one hand, the set of counterterms derived with a different method in [21] 
and, on the other hand, derive explicitly finite propagator and field equations. In Sec. IV we give the algorithms for 
solving at zero temperature and in Euclidean space the set of finite equations and the equations containing explicitly 
the counterterms. We demonstrate that their numerical solution coincide and investigate the efficiency of the two 
methods. We show that in the case of solving the equations containing the counterterms faster convergence of the 
iterative algorithm is obtained if the counterterms are rederived in order to allow them to evolve during the iterative 
procedure. Sec. V is devoted to conclusions. 

II. THE 2PI EFFECTIVE POTENTIAL AT TWO-LOOP FIELD-DEPENDENT LEVEL 

The 2PI effective potential of the real one-component <ft 4 model at the field-dependent two-loop truncation level 
can be given in the broken symmetry phase, characterized by the vacuum expectation value v, in the following form 



1 This distinction was introduced in [18] to differentiate the condition imposed one the 2PI vertex function obtained as field-derivatives of 
the inverse propagator which extremizes the 2PI effective action through the stationarity condition, from the condition which is imposed 
on the 2PI-rcsummed vertex function defined as field derivatives of the 2PI-resummed (1PI) effective action, and which differs from the 
2PI vertex function only due to the approximation made on the 2PI effective action. 

2 It was proven analytically in [30] that in Hartree approximation the phase transition cannot be of second order. 





FIG. 1. The counterterm diagrams corresponding to the terms of Vct[w, G] defined in (2). A wiggly line represents v, while a 
plain one corresponds to the full propagator G. Symmetry factors are not indicated. 




[6, 17, 21, 29]: 

V[v,G] = imV + A w 4 _ I J [ lnG -i (P) + jD - 1 (P)G(P)] +Uj p 

; / f G(K)G(P)G(K-P) + V ct [v,G}. (1) 
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The expression above is written in Minkowski space following the conventions of Ref. [31]: the tree-level propagator is 
D(P) = i/(P 2 — m 2 — Xv 2 /2) with P = (po, p) denoting the four-momentum and the metric is such that P 2 — Pq — p 2 . 
The counterterm functional is given by 3 



vy„, o, _ I Sm y + + ' fa + ^jc {P) + ^(jc m }\ 



(2) 



where the different terms are represented graphically in Fig. 1. As discussed in [17], the origin of two different mass 
counterterms and three coupling counterterms stems from the fact that within the 2PI formalism it is possible to define 
two two-point functions and three four-point functions, which do not necessarily coincide within a given truncation of 
the 2PI effective potential. In writing (1) it is implicitly assumed that, despite the possible different bare masses and 
counterterms their finite part agree, that is there is only one renormalized mass parameter m 2 , which in the symmetry 
breaking phase is usually negative, and only one renormalized coupling parameter A. 

The stationarity conditions SV[v, G]/5G = and SV[v,G]/Sv — give a self-consistent equation for the full two- 
point function and the equation for the vacuum expectation value, i.e. the field equation: 

iG^iP) = P 2 -£(P), (3) 

/ 1 1 \2 \ 

= v ( m 2 + 5m 2 l + -(\ + 5\4)v 2 + -{\ + 5\ 2 )T[G} + —S(0,G)\ , (4) 

where the self-energy S(P) is given by 

E(P) = to 2 + 5m 2 + i(A + 5\ 2 )v 2 + i(A + 5X a )T[G} + ^X 2 v 2 I(P, G). (5) 

Note that we call E(P) the self-energy, although its usual definition used in the 2PI formalism does not contain 
the mass parameter and the corresponding mass counterterm. The tadpole T[G], bubble I(P, G) and setting-sun at 
vanishing external momentum S^O, G) which appear above are defined as follows: 

T[G] = f G(K), (6a) 
Jk 



I(P,G) = -i [ G(K)G(P-K), (6b) 
Jk 



S(0,G) = -i G{K)G{Q)G{-K -Q). (6c) 
Jk Jq 

In [21] the divergences of these integrals were calculated by expanding the full propagator around the auxiliary 



3 Unfortunately in Ref. [21] Sm^ and <5m| were interchanged compared to the notation originally introduced in Ref. [17]. In this work we 
will stick to the notation of Ref. [21]. 
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propagator 



Go(P) = (7) 



with the mass scale Mo playing the role of the renormalization scale. The countcrterms absorbing these divergences 
were obtained in [21] by requiring the separate cancellation of the divergent coefficients of v , v 2 and of the environment 
(T and v) dependent finite function Tp[G], representing the finite part of the tadpole integral, both in the propagator 
and field equations. They are given in Appendix B in terms of some explicitly calculated integrals. 



III. DERIVATION OF EXPLICITLY FINITE EQUATIONS 

In this section we adapt the renormalization method used at finite temperature in [12] to the broken symmetry 
phase at zero temperature. We mention that the renormalization of scalar models displaying spontaneously broken 
symmetries was discussed generally in the context of the 2PI approximations in [16]. Here, we separate the field 
dependence from the divergent quantities in the same way as the temperature dependence was separated in [12] and 
give the renormalization conditions which, when imposed on the divergent quantities, lead on the one hand to exactly 
the countertcrms determined in [21] and summarized in (Bl), and on the other hand to explicitly finite equations, 
with no reference to any counterterms at all. We note that our countcrterms have exactly the same structure as those 
used in [29] and which were derived using the renormalization method of [17]. 4 In the next section we will check 
explicitly that the solution of the propagator and field equation obtained using counterterms agrees with the solution 
of the explicitly finite equations to be derived below. 



A. Renormalization of the propagator equation 



We begin by splitting the propagator into "vacuum" and "matter" parts, 

G(P) = G (v) (P) + G (m) (P), (8) 

where the vacuum part refers to a propagator defined in the symmetric phase (v = 0), while the matter propagator 
includes all explicit and implicit dependence on v. The truncation of the 2PI effective potential studied in this work 
leads to a momentum independent self-energy for v — which means that the vacuum propagator can be simply 
parametrized in terms of an effective positive mass M , as in (7), so that G^(P) = G a {P). Note that we need 
this propagator with positive mass squared M 2 because in the broken symmetry phase m 2 , the renormalized mass 
parameter of the Lagrangian could be negative. 

As a consequence of the splitting (8), the self-energy (5) is decomposed into three parts: 

E(P) = E (v) + E (0) (P) + E (r) (P), (9) 

where the explicit expression of the different parts are given below. The vacuum part is independent of the momentum, 
depends only on the vacuum propagator and has divergence degree 2. The last two pieces of the self-energy appearing 
in (9) are v dependent and are defined as follows. The part which does not result in any divergence when the full 
propagator is expanded around G^ is called regular part and is denoted by E^, while the remaining part, having 
divergence degree 0, is denoted by E^ ). The splitting of the field-dependent part is somewhat arbitrary. Using (8) in 
(5), the different parts of (9) are identified as 

EM = Ml + 5ml A + f GV{K), (10a) 

E(°)(P) = m 2 - Ml + Sml B + ^a£\o, P) + J^HK), (10b) 

E (r) (P) = -»AV / G {m) {K)G {v) {K + P)-i^- J G (m) {K)G {m) {K + P), (10c) 
Jk 2 J K 



4 Unfortunately the counterterms used in [29] appeared only on the corresponding poster. 
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where a£\o, P) is defined as 

A< v) (0,P) = A + <5A 2 -i\ 2 I G^{K)G {V) {K + P). (11) 



K 



We will see that with the above choice for the different pieces of the self-energy the expression of the countertcrms 
(Bl) can be obtained through simple renormalization conditions. The mass counterterm 8m\ was split into two parts: 
8m\ A is responsible for the renormalization of the vacuum part, while 8m\ B has to remove the divergence generated 
by the £(°) dependence of C?( m ). In order to see how this latter divergence proportional to £(°) emerges, we expand 
the propagator around the vacuum part: 

r(p] = * ( y , SW(P) + EW(P) 

1 > p2_ S (v) _ s(0) (P) - E« (P) P 2 - £< v ) (P) V P 2 - E< v ) 

= G (v) (P) - i (G (v) (P)) 2 E (0) (P) + G (r) (P). (12) 

The term proportional to E^ - 1 is C(l/P 4 ) (up to logs), therefore it gives divergent contribution after integrating over 
the momentum, while the regular part of the propagator goes with C(l/P 6 ) and leads to a finite contribution upon 
integration over the momentum. The sum of these two terms is the previously introduced matter part 

G (m) (P) = -«(g (v) (P)) 2 £(°>(P) + G (r) (P). (13) 
As already announced at the beginning of this subsection, the first renormalization condition is 

E< v ) = Mg , (14) 

which determines 8m\ A through the relation 



5mi A = -^Tf\ (1.1, 

where ' is defined in (B2). 



Now we turn to the renormalization of £(°). After plugging (13) into (10b), one obtains the following integral 
equation: 

E(o)(P)=m 2 -M 2 + ( 5< B + ^Ar ) (0,P)-,^l ^ ( G (v) (jr) ) a E (o) {K) + hi*** ^(K). (16) 
One can search for the solution of E^ ' in the following form: 

EW(P) = r« + y4 v) (0,P) + IrM J G^{K). (17) 

Using (17) in (16) one finds that Tm\ T^^O, P), and r' v ) fulfill the following Bethe-Salpeter-type equations: 

r« v) = ™ 2 - Mi + Sml B J k {g™[K))\ (18a) 

rW(o,p) = AW(o,p)-^ J (cM(i()) 2 rW(o,if), (18b) 

rM = A + S\ i±±*h r M J ( G W m ) 2 , (18c) 

in which all integrals are divergent. To render these divergent equations finite, one has to impose some renormalization 
conditions which will also determine the appropriate counterterms. Because of our introduction of the M scale and 
the splitting of the mass counterterm into two pieces, we need one more renormalization condition in addition to (14) 
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(ra . - ^ . - + (M . - «f - ^ (ra . - „>>f| / ( GW(J0 )», (1 „ 



to fix also 6m\ B . By using (10a) one sees that (m 2 — Mq) " 2 satisfies 

(jlvl ( 



which is the same equation fulfilled by r„\ up to counterterms. Since OH^/OMq = 1, a simple renormalization 

lt> = m 2 - Ml (20) 



condition on Fm is 



This fixes the 5m\ B counterterm: 



Sml B = (m 2 - M 2 )^!^ = -A±^£ (m2 _ M 2 )Tf , (21) 



where is defined in (B2). From (15) and (21) one obtains the complete 5m 2 . counterterms, as the sum of (5m 2 A 
and 5m 2 B , with the expression given in (Bla). Before moving to the renormalization of (18b) and (18c) we note that 
<5m| can be equivalently fixed to the same value by keeping the condition (20) but replacing (14) with the condition 

where "overall" refers to the environment independent part of the self-energy, which does not coincide with EW due 
to additional terms coming from and contains the usual overall divergence. One can identify this part of the 

self-energy as the one which remains after v and G ( T > are formally set to zero. Then, using the condition (20) in the 
right-hand side of (18a), with the notation introduced in (B2) one obtains in (17) E(°)(P)| = m 2 — Mq + 5m\ B + 

(A + 5A )(m 2 - M 2 )pj 0) /2, so that upon adding it to T,^ of (10a), the condition (22) becomes 

m 2 + 6ml + (pf + (m 2 - M 2 )Tf ) = m 2 . (23) 

Next, we focus on the renormalization of Ij an d r' v ^. As a renormalization condition we impose 

r< v » = A, (24) 
which using (18c) leads to (Bib) for 5A . Then, using (18b) and (11), one sees that the difference 



r 2 v) (o,p)-ri v) (o,o) = -u 2 



K 



G (v) (K + P)G (v) (K) - (g (v) (#)) ' 



A 2 4 V) (P) (25) 



is finite, where we defined Ip^ (P) as the finite part of the bubble integral in the zero-momentum subtraction scheme 

in which Ip\o) = 0. This shows that one only needs to renormalize the momentum independent 1^(0,0) which is 
achieved by imposing on it the renormalization condition 

4 V) (0,0) = A. (26) 

Then, the explicitly finite 1^(0, P) function is given by 

r 2 v) (0,P) = A + A 2 4 v) (P). (27) 

Using (27) in (18b) one obtains the equation determining <5A2, that is (Blc). 



In conclusion, the imposed renormalization conditions give explicitly finite expressions for the functions r„ , 
2^(0, P), and r( v ) appearing in (17), and provide also the following finite equation for E(°)(P): 

S (0) (P) = m 2 - Mq + y (A + A 2 4 v) (P)) + ^ J G^{K). (28) 
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FIG. 2. Diagrammatic renormalization of the « 2 -dependent part of the self-energy. See the text for details. The symmetry 
factors are not indicated and we used the shorthand Ao = A + 5Ao. The boxes denote the divergence of the encircled part of the 
graph, the plain line denotes the propagator G^ V '(P) = Go(P), while "F" denotes the finite part of the bubble integral. The 
index I in the sums goes over the possible number of bubbles (for the two classes of diagrams the lowest number of bubbles is 
zero or one, as indicated). 



Two interesting remarks are in order at this point. The first concerns the solution of (18b) which, as one can check 
iteratively, is easily expressed in terms of r' v ) as 

T^(0,P)=A 2 v) (0,P)- l -T^ J k A^(0,K)(g^(K)) 2 . (29) 

The second comment refers to obtaining a different expression for the 5X 2 counterterm. Using (24), (26), and (27) in 
(29) one obtains 

= SX 2 + iA(3A + SX 2 )T^ + ^A 3 ((Tf } ) 2 + T&) . (30) 

Note that (30) does not coincide with (Blc), but the solutions for 5X 2 do. Indeed, using the result of (Bib) for <5A in 
(Blc), the same expression as that coming from (30) is obtained. The Eq. (30) served in [21] as a consistency check 
of the renormalization procedure. 



We close this subsection with a diagrammatic illustration of the renormalization of the part proportional to v 2 in 
Eq. (10b) for £(°) , see Fig. 2. This is interesting because it shows the role the coupling counterterms play in the 
removal of both subdivcrgences and overall divergences. The first line in the figure shows the ^-dependent diagrams 
generated upon iterating Eq. (10b) after the Bogoliubov-Parasiuk-Hepp-Zimmermann (BPHZ) subtraction procedure 
was implemented (see [17] for details.) The counterterm <5A 2 displayed in Fig. 1 is decomposed in the sum of two 

terms SX 2 PUZ and AA 2 , and SX 2 PHZ is used to absorb the divergence of the bubble integral present in K^\q,P). 
The remaining finite part of the bubble integral is denoted by 'F'. In order to obtain the second line of the figure we 
localized all the subdivergences. Then, we used the renormalization condition (24). To obtain the third line of the 
figure, a cancellation occuring between the third and fourth term of the second line was exploited and the consequence 
of the renormalization condition (24), that is A — Ao = (AoA/2) j K G^\K), where Ao = A + 5Xq, was rewritten in a 
diagrammatic form. The remaining overall divergences of the diagrams in the third line are removed by AA 2 . Since 
the divergent integrals associated with these diagrams are and keeping in mind that one has to add also 

(5Af PHZ , from the condition of vanishing of the square bracket in the last line of Fig. 2 one obtains the expression of 
5X2 determined previously. 
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B. Renormalization of the field equation in the broken symmetry phase 

We use in the field equation (4) the splitting of the propagator introduced in the previous subsection in (8) and 
obtain 



Q = m 2 +5m 2 + ±±^ f G(v)( P )_A f f G<r\P + K)GW(P)GW{K) 
2 Jp 6 JpJk 

v 2 + \J K { ^{Q,K)G^\K) -i^r J J (3G {v) (P)+G (m) (P)^G {m) (P + K)G {m) (K). (31) 



A r OA 



6 

Using (13) and (17) in the second line, we obtain 



= m 2 + 5ml + / GW(P)-i- / / G^ V \P)G^(P + K)G^(K) 

2 J P 6 JpJ K 

J k ^\q,K)(g^{K)) 2 + V \^ + 1 J^\0,K)G^(K) 
-iy J J (3G (v) (P) + G (m) (,P))G (m) ( J P + ^)G (m) (^), (32) 
where we used the expression of given in (29) and defined 

if' = X + SX^-^z J^K 2 {Q,K){g^{K)) 2 Y^{0,K). (33) 



The integral in the expression above is divergent, and one imposes the following renormalization condition on T 



(v) 
4 



ri v) = A. (34) 

This leads to the equation for 6X4 given in (Ble). 

Since the last but one term of (32) is already finite in view of (26) and (27), and the last term is a finite integral, 
we are left with the divergences of the first line and the first term of the second line. In order to complete the 
renormalization one only has to find the quantity on which the renormalization condition determining 5ml can be 
imposed. It is straightforward to check that 



S 2 V[v,G v ] 



5v5v 



m 2 + 5m 2 + J G V = (P) - 1^ J ' J " G V=0 (P + K)G V = (P)G V=0 (K), (35) 



where G v on the left-hand side is the solution of the stationarity condition 5V[v,G v ]/5G v — and the subscript 
indicates explicitly that it depends on the vacuum expectation value v, so that the chain rule has to be applied 
when taking the derivatives with respect to v. With the splitting of the propagator G v =o into vacuum, matter, and 
"regular" parts, one can proceed exactly as at the beginning of this subsection, only that one has to omit everywhere 
the ^-dependent terms. Then, one has 

S2V ^ G ^ =m 2 + Sm l + ±+3l f G (v) (P) _ z A ! f f gW(P)GW(P + K)GW{K) 

2 Jp 6 JpJk 

J K ^\^K)(G^\K)) 2 + \ J t£\0,K)G<;%(K) 



5v5v 



v=0 



.A 2 
~ l ~6 



J p J (*Gi%{P) + G^l(P)) G<rt(P + K)G^l(K). (36) 



The last term above is finite and since T^\o,K) was already made finite by imposing the condition (26), the last 
but one term above is also finite. Therefore, one could in principle formulate our last renormalization condition as 

S ^SvSv^ lv-o = 777,2 ' then 5ml determined through this condition would differ by a finite term from 5ml given in 
(Bid). In order to obtain the desired expression for 5ml one choose the following condition 



5 2 V[v, G v 



5v5v 



= m 2 , (37) 

v— 0, overall 
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where the concrete expression on the left-hand side is given by the right-hand side of (36) without the last two terms, 
and the overall notation indicates that the divergence of this term is the usual overall divergence, that is the divergence 
which remains after all subdivergences are removed. 5 Then, it is easy to check that (37) gives indeed for <5toq the 
expression in (Bid). 

The two renormalization conditions (33) and (Bid) together with the expression of T^\o,K) given in (27) leads 
to the explicitly finite version of the field equation 

o = ™ 2 + ^ 2 + ^ J (l + XlP(P))G^(P) 

-i— [ [ G^(P + K)G {ra \P)G {ra \K) - i— [ [ G^(P + K)G^\P)G^\K). (38) 
2 JpJk 6 JpJ K 

The renormalization conditions imposed above give exactly the same expressions for these counterterms introduced 
in the functional (2) as those appearing in [21] (see also [Bl)]. This demonstrates the equivalence between the 
two renormalization procedures. We note that even though we formulated our set of renormalization conditions 
on the quantities T^ v \ 1^(0, P), corresponding to the three four-point functions V, V(0, P) = 5 2 T,(P)/Sv 2 , 

V = S 4 V[v, G v ]/5v 4 , introduced in [17], and also using the curvature S 2 V[v, G v ]/SvSv, they do not look as natural as 
those imposed in [17, 29] at T = T* and in the symmetric phase of the model (v = 0). In our case G (v) is not the full 
propagator at v — and in order to obtain a given expression of the counterterms which contains G ^ in place of the 
full propagator we need a specific way to formulate the renormalization condition in (22) and (37). One may even 
say that some of our conditions used for renormalization are not genuine renormalization conditions because they 
are not formulated on quantities accessible from the effective potential through functional derivatives and resemble 
more the way divergences are removed in a minimal subtraction scheme. As we have already formulated, we imposed 
these conditions in order to facilitate the comparison, to be done in the next section, between the solution of the 
finite equations and those using counterterms. On the other hand, as discussed in Ref. [30], finding a set of natural 
renormalization condition, prescribing e.g. the mass defined from the variational propagator and the curvature at the 
minimum of the effective potential 6 together with the value of the independent four-point functions is problematic 
when renormalizing the model at T — 0, that is in the broken symmetry phase, even at the Hartree level approximation 
of the 2PI effective potential. 



IV. NUMERICAL IMPLEMENTATION AND RESULTS 



Before discussing the algorithms and presenting the results we make some general statements about the numerical 
method used. 

Since we solve iteratively the set of coupled integral equations consisting of the self-consistent propagator equation 
and the field equation, one has to store and upgrade in each iteration step the value of some functions on a grid and 
to approximate using interpolation the value of these functions at the points required by the integration routine. For 
this purpose we use an equidistant grid for the modulus of the four-momentum, the one-dimensional Akima spline 
interpolation method and the numerical integration routines of the GNU Scientific Library [32] . When we solve the 
set of finite equations, one stores the regular and the matter part of the propagator [see (13)], while when solving the 
set of equations with the counterterms it is the full self-energy which is stored on the grid [see (9 and (10)]. In the first 
case the integrals are convergent in the UV, but nevertheless, for the functions which are stored on the grid we have 
to take into account the fact that there is a maximal value of the modulus of the four-momentum. The convolution 
of these functions should be calculated accordingly, using functions to restrict the momenta of the propagators, 
see (Al). For functions of the momentum with a known analytical expression there is no need to store them on the 
grid, and in consequence a simplification is encountered in the convolution involving such functions, see e.g. (A3). 
In contradistinction to this, when counterterms are used, all momentum-dependent functions are cut at the maximal 
value of the modulus of the four-momenta, that is at the physical cutoff A, irrespective of the fact that they are stored 
or not stored on the grid, in the latter case their expression is known analytically. We mention here that cutting the 
sum of momenta in the case of the bubble integral and as a consequence in all the double integrals encountered seems 
to be the cleanest way to proceed. 7 In fact, in the case of solving the set of finite equations we could afford to simplify 
the calculation, by cutting only the loop momenta in the integrals involving the propagator G^(P) = Gq(P) which 



5 The quantity on the left-hand side of (37) is obtained formally by setting v and G ' r ' to zero after (8), (13), and (16), (17) are used in 

the explicit expression of S ^^f^ ■ 

6 This two quantities do not necessarily coincide in a given truncation of the 2PI effective potential. 

7 See [33] for a discussion on the way such a regularization can be implemented at the level of the path-integral. 
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is not stored on the grid, because the integrals are all convergent in the UV, and in consequence the corrections are 
suppressed. We will sec indeed that the solution of the finite equations nicely agrees with that obtained by solving the 
equations with counterterms. Moreover, cutting only the loop-momenta, when solving the set of equation containing 
the counterterms, produces some small differences compared to the case when the sum of momenta is also cut, at least 
up to the largest cutoff we investigated, and in consequence is not a viable method for obtaining accurate numerical 
results. 



A. Algorithm for solving the finite propagator and field equations 



We will solve in Euclidean space the explicitly finite field and propagator equations obtained in the previous 
section. A Wick rotation is performed in every integral appearing in our equations: for example with the continuation 
ko — > ik 4 the tadpole integral reads J K G(k , k) — > i G(ifc 4 , k) = J feE A(k E ) where the Euclidean four-momenta is 
kE = (k, fc 4 ), such that k\ = k 2 + k\ and the Euclidean propagator is defined as A(k E ) = l/(k E + H(k E )). 

Using the expansion (12) of the propagator around the vacuum propagator A^I^e) — ^-/{k% + M 2 ) and the 
definition (13), the matter and regular parts of the Euclidean propagator reads 

A '°" <t£) = k% + AfJ + Z(<»(k E ) + EM(JL E ) " k% + Ml (39) 
A« (M=a <., (trf + _|^|_. (40) 

The two pieces of the self-energy which appear above are given by 

£(°>(M = m 2 -M 2 + A(l + A4 v) (M)y + ^ / AW(* B ), (41) 

\ 2 2 

E«(fc B ) = \ 2 v 2 I^ v \k E ) + -^-I (mm \k E ), (42) 
where the "matter-vacuum" and "matter-matter" bubble integrals are defined as 

I {mv) {k E ) = - I A^(pe)A^( P e + k E ), (43a) 
I {mm) (k E ) = - f A^(p E )A^(p E + k E ). (43b) 



The finite part of the vacuum bubble integral, I E \ appearing in (41), is defined in (B3). After analytical continuation 
to Euclidean space it has the following expression (recall that Go = G ^ ) : 



r(v) r * 1 



i ( ;>(k E ) 



16tt 2 



4M 2 



4M, 



(44) 



Note that a self-consistent equation for the regular part AW can be obtained in principle if one substitutes the 
expression of A^ m ^ from (39) into (40) because S' ' contains the integral of A^ and A^ appearing in the integrals 
of E« can be expressed in terms of A^ r ^ and Nevertheless, we store on the grid two quantities, A' m ) and A^\ 
which have to be solved simultaneously with the field equation (38), which in Euclidean space reads 

^ 9 A n A 



6 2 . lkE 



P \ 2 

/ (l + XlP(k E ))A^(k E ) + — (35 (mmv) + 5' (mmm) ), (45) 

Jk E V 1 6 

where the corresponding pieces of the setting-sun integral are defined as 

5 (mmv) = -/* f A^(p E )A^(k E )A^(p E + k E ), (46a) 

5 (mmm)^_/* f (p E ) A^ (k E ) A^ (p E + k E ) . (46b) 

Jpe Jkf- 
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Our iterative algorithm for solving the coupled set of equations for A^\ A^- m \ and v goes as follows. One starts at 
zeroth order by neglecting all the integrals in the self-energy and the field equation, including the finite part of the 
bubble 4 V) - By doing so, from (45), (39), and (40) one obtains: 



6m 2 (m) m 2 + ±v 2 -M 2 (r) (m 2 + \v 2 - M 2 ) 



A v &; (jfe| + Af^)(jfe|+m 2 + |w 2 )' ° (fc| + M 2 ) 2 (fc| + m 2 + |v 2 )' 

Then, starting from the /irsi order, the iteration is done in two steps. At a generic order n > 1 one first upgrades on 
the grid and A$ using (39) and (40), respectively, where the self-energy, which includes now Ip\ is calculated 
with the lower order quantities: v n -i, A^ m \, A^_i- As a second step of the nth iteration, one upgrades the vacuum 

expectation value by calculating v n from (45) with the integrals done with the already upgraded a!" 1 ' and An^ . This 
procedure is repeated until the iteration converges. The angular integration in the integrals appearing in (42) and 
(45) are performed in Appendix A as outlined at the beginning of this section. 

B. Algorithm for solving the propagator and field equations using counterterms 

If one decides to solve the propagator and field equations using counterterms, whose role in this equations is to 
cancel the divergent part of the integrals, then one needs a method to determine them. Such a method was developed 
in [21], and we review it here because it will be slightly modified for numerical reasons. In terms of the remaining 
finite part of the tadpole T F [G] and bubble integral I F (P, G) the propagator equation (3) reads 

iG- 1 (P) = P 2 -M 2 -^\ 2 v 2 I F (P,G), M 2 = m 2 + ^v 2 + ^T F [G}. (47) 
Using the expansion of G around Go 

G(P) = G (P) - iG 2 (P) (V 2 - M 2 + ^AV/ F (P, G)^j + ... , (48) 

the tadpole, bubble, and setting-sun integrals can be explicitly decomposed into divergent and finite parts (see Ref. [21] 
for details): 

7(P,G)=Tj 0) +I F (P,G), (49a) 
T[G] = rj 2) + (M 2 - M 2 )T^ 0) + ^AVTj J) + T F [G], (49b) 

5(0, G) = S (0) + 3Tj 0) T F [G] + 3(M 2 - M 2 ) ((rj 0) ) 2 + rj J) ) + |aV (t^T^ + Tf 2) ) + S F (0, G), (49c) 

where the divergent integrals T d (2) , Tj 0) , T { J ] , T^' 2) and 5 (0) arc defined and calculated in Appendix B. Plugging 
the decomposed expressions of the integrals (49) into the propagator and field equations and subtracting from them 
the corresponding explicitly finite equations, written in terms of the renormalized parameters m 2 ,A and the finite 
part of the integrals, one obtains relations between the counterterms and the divergent integrals. Requiring that the 
coefficient of the v 2 , v°, and Tp[G] vanish independently, one obtains the counterterms determined in [21]. These 
were also obtained from appropriate renormalization conditions in Sec. Ill and are summarized in Appendix B. 8 

In Euclidean space these counterterms are functions of the four-dimensional rotational invariant cutoff A appearing 
through the integrals explicitly calculated in (B5). They can be used to solve iteratively the set of Euclidean equations 

S(p) = m 2 + Sm 2 + 1(A + 5X 2 )v 2 + 1(A + <5A )T(£) + ^X 2 v 2 I(p, S), (50a) 

= m 2 + Sm 2 + 1(A + S\ 4 )v 2 + h\ + <5A 2 )T(£) + ^-5(0, S). (50b) 
6 2 6 

However, since the counterterms determined from (Bl) are designed to cancel in these equations the divergences of 
integrals produced by the full (i.e. iteratively converged) propagator, they cancel the complete cutoff dependence of 



Actually, applying the outlined procedure to the field equation one obtains (30) from the cancellation of the coefficient of v 2 , but as 
stated below (30), its solution for <5A2 coincides with the one obtained from the propagator equation, as it should. 
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T[£], I(P, E) and 5(0, E) only when the solution of this set of equations has converged. This is not a problem in itself, 
since, as it will be demonstrated, the iterative procedure converges, but the convergence of the solution of the set of 
equations is rather slow. It is possible to improve the iterative procedure if, instead of using the counterterms given in 
(Bl), we rederive them at each order of the iteration using the procedure outlined above. These counterterms, which 
come from the requirement to have finite equations at each order of the iteration, will evolve during the process of 
iteration toward the value of the counterterms obtained from (Bl) which will be reached when the iteration converged. 
Since, as we will see, the new counterterms are now environment dependent, that is depend on v and IV [E], this 
convergence can only happen if v and E also converge to their appropriate values, and has to be checked numerically. 

In what follows we present the improved iterative procedure which uses evolving counterterms. At each order the 
iteration starts by the upgrade of the self-energy, followed by the upgrade of the vacuum expectation value determined 
from the field equation. 

At zeroth order of the iteration the field and the self-energy are the tree-level ones: 

«2 = -^T> E (p)=™ 2 + ^ 2 . (51) 

At this order, since there are no quantum fluctuations at all, there are no divergences to cancel and in consequence 
all the counterterms are zero, that is SrriQ = 5m\ = <5A ,o = S\ 2 .o = <5A 4 , = 0, where the counterterms also carry 
the index of the iteration number because all of them will change during the iteration process. 

The general formulas for the nth order counterterms to be given below encompass all orders with the following 
prescription: at n = one starts from the tree-level values, and formally T(E„) = /(£„) = for n E { — 1,-2} 
and V-i = vq. Suppose that the (n-\)th order of the iteration is done, i.e. v„_ 1 and E„_i(p) together with the 
corresponding counterterms are already known, then at nth order we start by upgrading the self-energy using [c.f. 
(50a)]: 

£„(p) = m 2 + Sml n + \vl_, + ^ 2 _ 2 + ^±|^T(E n _ 1 ) + J(p, £„_>). (52) 

The finite part of the tadpole and bubble are defined through 

T(S n _ 1 ) = if + (m 2 + \vl_ 2 + ^T F (E„_ 2 ) - M 2 )rf + ^f^T« + 2>(E„-i), (53) 

I(p, E„_ x ) - lj 0) + I F (p, E n _i). (54) 

The renormalization procedure requires the independent vanishing of the overall divergence, field-dependent subdi- 
vergence and TV-dependent subdivergence: 



= 8m 2 



n + ^±^(TP + (^-M^), (55a) 



\ \ 2 

= 5\ 2 . n + -lf\\ + SX , n ) + -Tj 7) (A + <5A ,„)(1 - S n ,i) + A 2 rj 0) ^p±, (55b) 

A 1 V n-2 

= (A + ( 5A ,„)^ri 0) T F (E„_ 2 ) + 5Ao, n T F (£„_i). (55c) 

The counterterms (5m 2 , „, 6X2^ and (5Ao,n are obtained from (55), so that E n (p) is finite and can be calculated. 
Note that for n = 1 the second term in the equation for <5A2.„ vanishes as it should because in this case there is no 

divergence due to the bubble, that is there is no T ( d I] in T(E ), whose expression is 

T(E ) = Tf + (m 2 + - M 2 )rf + 7>(E ). (56) 

Note also that for n = 1 one has <5Ao,i = 0, since Tp(E_i) = by convention. 

We continue with the determination of u 2 from the field equation, which at nth order of the iteration reads 

= m 2 + Sml n + ±vl + ^vl^ + A±^i T(S „) + * S (0, E„), (57) 
where T(E„) is given by (53) with n — 1 — > n and the setting sun is 

5(0, E„) = 5 (0) + 3Tf T F [E„] + s(m 2 + \v\_ x + ^T F (E„_!) - M 2 ) ((rf) 2 + T^) 
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+ ^A\ 2 l _ 1 (Tj J) Tf +Tf 2) ) +5 F (0,E„). (58) 
Applying the same requirements as for the self-energy, one obtains: 

= K„ + A±^i( T f + (m > _ ^ )T (0)) + * ( m » - M^(T^f + T«) + £s (0), (59a) 

= ^A 4 ,„ + ^(A + SX 2 , n + A 2 Tf ) (if + XTP) + ^ (l^ + A7f 2 >) , (59b) 

= a 2 ,„ + A 2 Tf + [£(A + *A 2 ,„)Tf + £ ((Tf ) 2 + if)] . (59c) 

Note that a new counterterm SX 2 has been introduced instead of SX2, since in the nth iteration SX2,n does not eliminate 
properly the appropriate subdivergences of the field equation. Extracting 5m 2 , n , SX4 yrl and 6X2^, from (59) one can 

obtain v 2 from (57). The counterterm 8X 2 , n will be equal to <5A2, n determined from the propagator equation only 
when the iteration converged, so that Tp(S„_i)/Tp(S„) — > 1. In this limit and when v n /v n -\ — > 1 the counterterms 
will converge to the values obtained from (Bl). This happens only asymptotically and in practice the iteration stops 
when the relative change in the values of the field and the self-energy are smaller than some given value dictated by 
our requirement of accuracy. The same iterative procedure is used when the counterterms are not evolved during the 
iteration, but used as determined from (Bl): one first upgrades the self-energy and then the field equation and checks 
for their convergence. 



C. Numerical results 



In what follows we present the results on the numerical solution of the propagator and field equations in Euclidean 
space. When solving the expli citly finite set of equations the maximal value of the modulus of the four-momentum 
stored on the grid is Lj y / 2\m 2 \ — 500, while in the case of solving equations with counterterms the highest value of the 
modulus stored, that is the four-dimensional physical cutoff, is A/ yj2\m 2 \ = 200. The value of the mass scale is chosen 
to be M / ^/2\m 2 \ = 2.1 and the value of the renormalized mass parameter squared is m 2 = —0.5, in some arbitrary 
physical units. Actually, to have an idea on the physical scale, one can require to have v/Y>^(p = 0) ~ 0.66, which is 
the value of the ratio f v /m v of the pion decay constant to the pion mass for f v — 93 MeV and — 140 MeV, and 
£2 (p = 0) = m n . Then, from the left panel of Fig. 3 one sees that the first condition can be meet by choosing A ~ 8, 
for which value the second condition means y / 2\m 2 \ ~ 158 MeV. In units for which m 2 = —0.5 the step size on the 
grid in general was 0.01. We have checked that by halving the step size the change in the results is comparable with 
or smaller than the precision required by our convergence criterion. The iteration was stopped when the change in 
the vacuum expectation value was smaller than 10~ 6 . The precision of the numerical integration routines used was 
higher than that, since we required relative error bounds between 10~ 6 and 10 -8 . 

The solution of the explicitly finite equations can be seen in Fig. 3. In the right panel one can see the converged 
solution of the field equation obtained for various values of the coupling as function of L, the maximal value of the 
modulus stored on t he grid . As explained in the figure caption, v was shifted appropriately in order to make all four 
curves meet at L/y/2\m 2 \ — 500. In all cases a plateau can be observed for increasing values of L. However, with 
increasing value of the coupling A the plateau starts at a higher L. With the same convergence criterion the number 
of iterations increases with increasing A and for A = 30 one needs twice as many iterations until convergence occurs 
compared to the A = 3 case. Interestingly, at a given A, the number of iterations does not depend on L. 

Next, we would like to quantify the extent of achievable improvement given by the use of the iteratively evolved 
coupling and mass counterterms, as compared to the solution of the equations which use the counterterms originally 
derived in [21]. In these two cases we compare in Fig. 4 the number of iterations needed for the convergence of v at four 
different coupling constants. We use a reasonably large cutoff A/ yj2\m 2 \ = 200, around which the convergent results 
are expected to show cutoff independence within some desired accuracy. One observes the advantage of evolving 
the counterterms even for lower values of the renormalized coupling. For larger values (A > 10) the improvement 
is significant, as one can see from the large iteration number needed by the algorithm using fixed counterterms. 
Therefore, it is expected in general that solving self-consistent equations iteratively in a fully nonperturbative regime 
with fixed counterterms is not efficient from a numerical point of view. For large couplings, A > 10, the number of 
iterations needed in case of using fixed counterterms are at least a factor of 2 larger than the number of iterations 
needed for the explicitly finite equations to converge. 

In Fig. 5 we compare for A = 5 the ways how a cutoff independent solution is reached using fixed and evolving 
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FIG. 3. The solution of the explicitly finite propagator and field equations 
of the coupling obtained for L 



Left panel: v/M and M = (p = 0) as function 
The converged value of v as function of the maximal value of the 

An appropriate constant was 



/2\m 2 \ = 500. Right panel: 
modulus of the momentum stored on the grid for different values of the coupling constant A. 
subtracted from the values of v obtained for A = 3, 10, and 20 such that the resulting values coincide at L/y / 2\m 2 \ = 500 with 
the value of v obtained for A = 30. 



1 

0.996 
0.992 
0.988 
0.984 
0.98 

0.55 
0.54 
0.53 
0.52 
0.51 
0.5 



evolving cts 
fixed cts 



evolving cts 
fixed cts 



X=3 



X=10 



0.775 
0.77 

0.765 
0.76 

0.755 
0.75 

0.745 

0.38 
0.36 
0.34 
0.32 



evolving cts 
fixed cts 



1 



k=5 



4 5 



evolving cts 
fixed cts 



V* /. 20 





8 10 



10 15 



20 



FIG. 4. A comparison between the convergence of the itera tive al gorithms using fixed and evolving counterterms for a fix value 
of the physical cutoff: A/^/2\m 2 \ = 200. The value of v/y/2 \m?\ is shown as function of the number of iterations for different 
values of the coupling. 



counterterms. The left panel shows that with increasing values of A a plateau in the v(A) functions is approached 
rather slowly when the fixed counterterms of [21] are used, despite the fact that the coupling is not large. Looking 
at the plot with values of A spanning a large range it seams that after several iterations a plateau is finally reached. 
However, the structure in v(A), better seen in the inset which uses a smaller A range, shows that the convergence is 
not uniform: more iterations are needed for larger cutoffs to get close to the true solution, that is the solution of the 
explicitly finite equations (horizontal line). This nonuniformity of the convergence in A is more pronounced for larger 
values of the coupling constant. Contrary to this, the right panel of Fig. 5 shows a uniform convergence in A when 
the iterative algorithm uses evolving counterterms. Since by the very construction of the evolving counterterms a 
plateau is observed already for moderate values of the cutoff at every order of the iteration, one sees again that with 
this improved approach the final, convergent result is obtained in a more efficient way. 
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FIG. 5. A comparison between the convergence rate of the iterative algorithms using fixed (left panel) and evolving counterterms 
(right panel) as reflected in the change of the vacuum expectation value at different values of the physical cutoff A. The horizontal 
line represents the solution of the explicitly finite system of equations obtained with Lj ' \j2\m?\ — 500, where L is the maximal 
value of the modulus of the momentum stored on the grid. The insets show the convergence at large iteration numbers. 




FIG. 6. The convergence of the full self-energy obtained for A = 5 with evolving counterterms for A/ y/2\m 2 \ = 200. The inset 
shows the behavior of the self-energy at large momenta. 



The convergence of the self-energy can be seen in Fig. 6. The function converges smoothly, i.e. for all momentum 
values the self-energy reaches its converged value almost simultaneously. Note however that the convergence for 
Pe — is slightly slower as compared to larger values of the momenta. 



V. CONCLUSIONS 



We presented an approach which can be used to obtain a very accurate numerical solution of the self-consistent 
propagator equation coupled to the field equation, both derived from the two-loop level approximation of the 2PI 
effective action in the <p 4 model. We showed that the renormalization method developed in [21] which removes the 
divergences using minimal subtraction is equivalent to imposing nontrivial renormalization conditions on the self- 
energy, the curvature of the effective potential and on kernels of Bethe-Salpeter equations related to these quantities. 
The use of renormalization conditions also allowed us to construct explicitly finite propagator and field equations. 

We compared the convergence of the iterative method applied to solve, on the one hand, the explicitly finite 
equations and, on the other hand, the equations which contains the explicitly calculated counterterms of [21]. The 
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very same results were obtained numerically which indicates the correctness of our analytic study on renormalization. 
It turned out that, especially for larger coupling constants (A > 10), the convergence rate of the algorithm which uses 
the counterterms is quite poor. This should not come as a surprise given that the counterterms were determined based 
on the asymptotic behavior of the propagator. Within an iterative procedure the asymptotics develops progressively 
and is reached only when the solution converged. Only then the divergences present in the equations match the 
counterterms constructed to cancel them, at intermediate steps of the iteration process the use of fixed counterterms 
results in an oversubtraction. In order to cure this, we managed to develop a different algorithm which rederives the 
counterterms at every step of the iterative procedure from the condition that the counterterms cancel the divergent 
part of the integrals at every order of the iteration. This method which uses evolving counterterms turned out to be 
very effective in obtaining the converged solution within a few iterative steps also for larger couplings. 

The presented approach and some of the methods applied in our present work may be directly extended, first to 
a zero temperature numerical solution of the model in Minkowski space, and then to a finite temperature solution 
of the model, where it would be interesting to investigate the order of the phase transition. We expect that the 
solution which can be obtained with the method described here could serve as a good benchmark for other, even more 
effective methods used in finite temperature studies (see e.g. [28, 29]). Some of these investigations will be presented 
in forthcoming publications. We expect that our methods presented here will represent a good basis for obtaining 
highly accurate and/or numerically controlled solutions of the 2PI approximation in more complicated theories as 
well. We plan to carry out the study of the phase transition in the O(N) symmetric model using large- N expansion. 
This is motivated by the fact that properly renormalized complete solutions of the equations derived in the O(N) 
model at next-to-leading order of the 2PI-1/7V expansion are missing from the literature. The renormalization and the 
thermodynamics of the O(N) model was investigated recently at the next-to leading order of the 1PI-1/JV expansion 
in [34]. 
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Appendix A: Integrals appearing in the finite Euclidean equations 



Since the functions appearing in the quantum correction integrals depend only on the modulus of the four- 
momentum, certain parts of the angular integral and in some cases the entire angular integral can be performed 
analytically. The angular integration for the integrals involving A( r ) is trivial, while the integral over the modulus of 
the four-momentum can be performed only numerically. As explained at the beginning of Sec. IV for the convolution 
of two functions stored on a grid one uses 

C L [f,g](k E ) = J ^f(\ PE \)g(\ PE + k E \)9(L-\k E \)0(L-\k E +p E \) 

I r L /■min(p+|k E |,.L) , - 

= FT7t/ d PPfiP) dqqg(q)jAp 2 q 2 - (p 2 + q 2 - k%) , (Al) 

O 71 " fc B JO J\p-\k B \\ 

where 9 is the step function, and for the modulus we introduced p — \pe\, Q = \&e +Pe\, L is the maximal value of 
the modulus stored on the grid. The q integral is a remnant of the angular integration in a four-dimensional spherical 
coordinate system and is obtained with the change of variable cos6*i = (q 2 — k 2 — p 2 )/(2kp). When the function g is 
A( v \ which is not stored on the grid, that is there is no need for the second 9 function, the entire angular integration 
can be done analytically, and with the change of variable t = tan(#i/2) one has 



/' 

Jo 



^ p 2 E + k 2 E + 2^\\k E \ C o & 9 1+ M 2 - 4lM ( 4 + Pl + Ml ~ v(4+,| + Mo 2 ) 2 -4^l) • (A2) 
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Using (A2) in the first term of (42) and for the setting-sun contribution to the field equation (45) which contains one 
vacuum propagator one obtains 

l(™\k E ) = -jg^i £ d PP (k 2 E +p 2 + Ml - ^{k 2 E +p 2 + M 2 f -4k 2 E p 2 ^j A^(p), (A3) 

g (nm,v) = __±_^ dkkA^ a \k) J dpp (k 2 + p 2 + M 2 -^J (k 2 + p 2 + M 2 ) 2 - 4fc 2 p 2 ) A< m )(p). (A4) 

For the second term of (42) and the setting-sun contribution to (45) which contains only matter propagators one uses 
(Al) and the results of the angular integrations are 

1 r L rmin(p+\k B \,L) , 

I^\k E ) = T / dppA^{p) dqqJAp 2 q 2 -(p 2 + q 2 -k 2 E ) 2 A^(q), (A5) 

i pL pL p min (p-\-k,L) 

S* (mmm) = ;r / dkkA^(k) dppA (m \p) dqq^4p 2 q 2 - (p 2 + q 2 - k 2 ) 2 A^ m \q). (A6) 



Appendix B: Cutoff dependence of divergent integrals 

For reader's convenience we summarize below the equations for the counterterms obtained in [21] and rederived in 
Sec. Ill by imposing appropriate renormalization conditions: 

=Sm 2 2 + 1(A + 5Xo) [TP + (m 2 - M 2 )T { p] , (Bla) 

0=SX + ^X(X + SX )T^\ (Bib) 

=5X 2 + Ia(A + <5A ) (rf + XTP) + A 2 Tf , (Blc) 

=Sm 2 + i(A + SX 2 ) [Tf + (m 2 M 2 )Tf ] + ^(m 2 - M 2 ) ((if ) 2 + if) + 1a 2 S (0), (Bid) 

=SX 4 + h(X + 5X 2 + A 2 lf ) (if + Alf ) + ^A 3 (if + Alf 2 >) . (Ble) 

The first three and the last two expressions come from the renormalization of the propagator equation and field 
equation, respectively. The divergent quantities if, if, and 5o(0) arc the vacuum parts of the tadpole integral 
(6a), the bubble integral (6b) at vanishing external momentum and the setting-sun integral (6c), all calculated with 
the propagator Go given in (7): 



T (2) 



J p G (P), (B2a) 

TT = % = -ij p Olin (B2b) 

S (0) = -i[ f G (K)G Q (Q)Go(K + Q). (B2c) 

JK JQ 

The remaining two divergent integrals if and if 2 "* are defined in terms of the finite part of the bubble integral 
(6b) calculated with the propagator Go, that is 

I , F {K) = ~ij (Go(P)Go(P + K)- Gl(P)). (B3) 



and read 



if = -if G 2 (K)IoMK) = -if G 2 (K)I Q (K) (if ) 2 - (B4a) 

JK JK 

T { J' 2) = -i[ G 2 (K)I 2 AK) = -i[ G 2 (K)I 2 (K) 2lf if - (if ) 3 . (B4b) 

JK JK 
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With a four-dimensional rotational invariant cutoff regularization in Euclidean space, the divergent integrals are 
computed as follows 



T (2) _ 



A 2 



16tt 2 



Ml , / A 2 



Ml 



T (o) = _j_ 
d 16tt 2 



A 2 



A 2 + Ml 



So(0) 



dk 



A; ;i 



T, 



-A 



T U,2) _ 



dk- 

A 



(fc 2 + Af 2 ) 
k 3 



c a [aw,a (v) ]W, 



87r 2 y "'"(fc 2 + M 2 ) 2 



dfc- 



C A [A^\A^](k)- (T t 
fc 3 - -> 



(0)x2 



(cf A [AW, A M](fc)) - 2T { d I] Tf ] - (T^) 3 , 



(B5a) 
(B5b) 
(B5c) 
(B5d) 



8^ 2 7 ~'"(fc 2 + M 2 ) 2 

where for C A [A( V \ A< v )](fc) one uses (Al) with A( v '(fc) = l/(fc 2 + Mq) and with L replaced by the physical cutoff A. 
The counterterms determined from (Bl) using the expressions in (B5) cancel the divergences of tadpole, bubble and 
setting-sun integrals, which in Euclidean space read as 

T[A] = ^L^ dkk 3 A(k), I(k,A) = -C A [A,A](k), S(0) = -^^ dkk 3 A(k)C A [A,A}(k). (B6) 
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